class: inverse,left, middle background-image: url(data:image/png;base64,#background.png) background-size: cover <img src="data:image/png;base64,#LOGO_DIPLOMADO.png" width="500px"/> ##Módulo 2: Estadística espacial y geoestadística ### Interpolación por Kriging Ordinario Javiera Aguayo T.<br> javiera.aguayo@pucv.cl<br> .large[<b><a href="https://www.pucv.cl/uuaa/site/edic/base/port/labgrs.html">LabGRS</a> | Octubre 2023</b>] <br> --- class: center,middle background-image: url(data:image/png;base64,#labgrs_logo.png) background-size: 35% --- ##Contenidos .pull-left[ 1) Etapa 1 Modelado Geoestadístico: AEDE 2) Etapa 2 Modelado Geoestadístico: Análisis Estructural - Variograma - Parámetros de un Variograma - Estimación del Variograma Experimental - Ajuste del Variograma Experimental 3) Etapa 3 Modelado Geoestadístico: Predicción - Interpolaciones por Kriging - Interpolación por Kriging Ordinario (KO) - Varianza de la interpolación por Kriging Ordinario ] .pull-right[ <right><img src="data:image/png;base64,#https://ucsbcarpentry.github.io/CustomDC-R/fig/r_rollercoaster.png" width="500px"/></right> ] --- ##Etapa 1 Modelado Geoestadístico: AEDE ###Revisión del supuesto de normalidad .pull-left[ **- Histogrma:** <center><img src="data:image/png;base64,#histograma_AEDE.png" width="550px"/></center> ] .pull-right[ **Coeficiente de Sesgo =** 0.045 **Kurtosis =** 2.21 **Test de Normalidad:** - Jarque - bera **->** p-value = 0.1421 - Kolmogorov - smirnov **->** p-value = 0.0278 - Shapiro - wilk **->** p-value = 0.0115 ] --- ##Etapa 1 Modelado Geoestadístico: AEDE ###Revisión de Estacionariedad y autocorrelación espacial .pull-left[ **Cálculo de Moran:** p-value = 0.493794 **Correlograma:** <center><img src="data:image/png;base64,#correlograma.png" width="400px"/></center> ] .pull-right[ **Mapa variográfico:** <center><img src="data:image/png;base64,#mapa_variografico.png" width="550px"/></center> ] --- #### Lectura y visualización de los datos ```r datos_temp_cauquenes <- read.csv("temp_promedio_cauquenes.csv", sep = ",") ptos_temp_cauquenes <- vect(datos_temp_cauquenes, c("X", "Y"), crs = "epsg:32719") # EPSG poli <- vect("Provincia_Cauquenes.shp") ``` <img src="data:image/png;base64,#DIPGEOPR_02_6_files/figure-html/unnamed-chunk-3-1.png" width="100%" /> --- ## Etapa 2 Modelado Geoestadístico: Análisis Estructural ### Variograma Es una función matemática que muestra como cambia la semivarianza entre dos localizaciones a distintas distancias, y permite analizar la estructura y la dependencia espacial o autocorrelación de los datos muestreados. -- ####**¿Qué es la semivarianza?** .pull-left[ Se define como la mitad del valor esperado de la diferencia al cuadrado entre los valores de una variable en dos ubicaciones geográficas distintas. _Fórmula Estimador Empirico del Semivariograma:_ _(Matheron, 1963)_ <center><img src="data:image/png;base64,#Formula_semivarianza.png" width="500px"/></center> ] .pull-rigth[ <center><img src="data:image/png;base64,#Mapa_con_distancias.png" width="550px"/></center> ] --- ## Etapa 2 Modelado Geoestadístico: Análisis Estructural ###Parámetros de un Variograma .pull-left[ **- sill:** semivarianza máxima de la variable de interés. **- nugget:** indica que hay una variación a distancias cortas (habitualmente asociada a un error de medición.) **- range:** distancia hasta la cual existe autocorrelación espacial (semivarianza es máxima y no sigue aumentando). Para definir un modelo de variograma, se utilizará el **paquete gstat** y en particular la función `vgm()` (Variogram Model). ] .pull-rigth[ <center><img src="data:image/png;base64,#semivariograma.png" width="500px"/></center> ] --- ## Etapa 2 Modelado Geoestadístico: Análisis Estructural ###Estimación del Variograma Experimental ```r ## Calculo de semivarianza para todos los puntos v <- gstat(formula = t_mean1~1, locations = ~X+Y, data = datos_temp_cauquenes) cloud_v <- variogram(v, cloud = T) plot(cloud_v) ``` <img src="data:image/png;base64,#DIPGEOPR_02_6_files/figure-html/unnamed-chunk-4-1.png" width="100%" /> --- ## Etapa 2 Modelado Geoestadistico: Análisis Estructural ```r ## Variograma experimental v_exp <- variogram(v, boundaries = c(1000, 3000, 7000, 15000, 20000, 30000, 60000, 100000)) # Intervalos de distancia plot(v_exp, plot.numbers = TRUE) ``` <img src="data:image/png;base64,#DIPGEOPR_02_6_files/figure-html/unnamed-chunk-5-1.png" width="100%" /> --- ## Etapa 2 Modelado Geoestadístico: Análisis Estructural ### Ajuste del Variograma Experimental Una vez listo el variograma experimental, este se debe ajustar en base a un modelo seleccionado según la interpretación del usuario desarrollador del variograma experimental. -- Algunos de los modelos más utilizados son: **- Modelo Esférico (Sph):** es uno de los más populares entre los modelos de semivariograma. Tiene dos características principales: un comportamiento lineal cerca del origen y el hecho de que en terminos de distancia se alcanza un rango en que ya no hay autocorrelación espacial, donde el semivariograma encuentra la meseta y después de esta se mantiene llano. -- **- Modelo Exponencial (Exp):** este modelo se parace mucho al modelo esférico, pero se diferencia en que para el mismo rango y meseta de un modelo esférico, el modelo exponencial alcanza el rango más rápidamente, es decir, a menor distancia que el modelo esférico. -- **- Modelo Gausiano (Gau):** tiene un comportamiento cuadrático cerca del origen y produce una autocorrelación espacial de corto rango. --- ## Etapa 2 Modelado Geoestadístico: Análisis Estructural ```r ## Modelos disponibles para el ajuste del variogrma experimental show.vgms() ``` <img src="data:image/png;base64,#DIPGEOPR_02_6_files/figure-html/unnamed-chunk-6-1.png" width="100%" /> --- ## Etapa 2 Modelado Geoestadístico: Análisis Estructural ```r ## Ajuste del Variograma experimental parametros_variograma <- vgm(nugget = 10, # Error aletorio asociado al muestreo psill = 70, # Semivarianza máxima de la variable range = 60000,# Distancia hasta la cual existe correlación espacial model = "Exp")# Modelo de Variograma - Exp (Exponencial) v_ajustado <- fit.variogram(v_exp, parametros_variograma, fit.method = 6) #Ajustes por minimos cuadrados class(v_ajustado) ``` ``` ## [1] "variogramModel" "data.frame" ``` --- ## Etapa 2 Modelado Geoestadístico: Análisis Estructural ```r plot(v_ajustado, 80000) ``` <img src="data:image/png;base64,#DIPGEOPR_02_6_files/figure-html/unnamed-chunk-8-1.png" width="100%" /> --- ##Etapa 3 Modelado Geoestadístico: Predicción ###Interpolaciones por Kriging El Kriging es uno de los varios métodos utilizados para realizar interpolaciones, la cuales constan de poder realizar estimaciones de valores no muestreados en base a un conjunto limitado de datos muestreados. -- ####Características principales: -- - A diferencia de los métodos deterministicos, **_el kriging utiliza la correlación espacial entre los puntos muestreados para interpolar los valores en el espacio._** -- - Los pesos asignados por el kriging, se calculan considerando, que los puntos cercanos a la ubicación de interés reciben más peso que los más lejanos. -- - Las ponderaciones consideran la disposición espacial de todos los datos muestreados, por lo que **_las agrupaciones en áreas sobremuestreadas se ponderan menos (menor peso)_**, porque contienen menos información que los datos individuales. -- - Este método es considerado el mejor **_estimador lineal no sesgado._** --- ##Etapa 3 Modelado Geoestadístico: Predicción ### Interpolación por Kriging Ordinario (KO) Este tipo de kriging se basa en los siguientes supuestos: -- - Se asume que el promedio de los datos muestreados es constante. -- - Se asume que existe dependencia espacial o autocorrelación espacial, es decir existe un varigrama, aunque este es desconocido. -- Para calcular el Kriging Ordinario, se utilizan los valores del variograma ajustado para obtener predicciones en ubicaciones no muestreadas. Para lo cual utilizaremos la función `gstat()` del paquete **gstat**. --- ##Etapa 3 Modelado Geoestadístico: Predicción ```r ## Creación de una grilla sin información raster_vacio <- rast(extent = ext(poli), #Extensión del poligono de entrada resolution = 500, #Resolución espacial crs = "epsg:32719") #Sistema de proyección raster_poli <- rasterize(poli, raster_vacio) plot(raster_poli, legend = F) ``` <img src="data:image/png;base64,#DIPGEOPR_02_6_files/figure-html/unnamed-chunk-9-1.png" width="100%" /> --- ##Etapa 3 Modelado Geoestadístico: Predicción ```r ## Interpolación por Kriging Ordinario Kriging_ordinario <- gstat(formula = t_mean1~1,#Formula kriging ordinario locations = ~X+Y, #Columnas X e Y data = datos_temp_cauquenes, model = v_ajustado) interpolacion_ko <- interpolate(raster_poli, Kriging_ordinario, xyNames = c("X", "Y")) %>% mask(raster_poli) ``` ``` ## [using ordinary kriging] ## [using ordinary kriging] ``` ```r names(interpolacion_ko) <- c('Interpolación', 'Varianza') ``` --- ##Etapa 3 Modelado Geoestadístico: Predicción <img src="data:image/png;base64,#DIPGEOPR_02_6_files/figure-html/unnamed-chunk-11-1.png" width="100%" /> --- ##Etapa 3 Modelado Geoestadístico: Predicción ### Varianza de la interpolación por Kriging Ordinario La función `gstat()`entrega dos resultados, el primero es la predicción y el segundo es la varianza de la interpolación por kriging ordinario, la cual nos permite cuantificar la incertidumbre de las predicciones. -- Este mapa de varianza será el primer indicador que utilizaremos para evaluar la calidada de la interpolación realizada. -- ```r #Varianza de la interpolación por kriging plot(interpolacion_ko[[2]], main = "Mapa de Varianza", col = rev(heat.colors(256)), legend=TRUE, zlim=c(0,60)) lines(poli) ``` --- ##Etapa 3 Modelado Geoestadistico: Predicción ### Mapa de Varianza <img src="data:image/png;base64,#DIPGEOPR_02_6_files/figure-html/unnamed-chunk-14-1.png" width="100%" /> --- ### Bibliografía 1963.G. Matheron. Principles of Geostatistics. Economic Geology, Vol.58. http://cg.ensmp.fr/bibliotheque/public/MATHERON_Publication_02396.pdf 2003.Rubén Fernandez Casal. Geoestadistica Espacio-temporal: Modelos flexibles de variogramas anisotrópicos no separables. https://rubenfcasal.github.io/files/Geoestadistica_espacio-temporal.pdf 2014.Edzer J. Pebesma. gstat user's manual. http://www.gstat.org/gstat.pdf 2022.Edzer J. Pebesma.Package ‘gstat’. https://cran.r-project.org/web/packages/gstat/gstat.pdf 2022.Rubén Fernández Casal y Tomás Cotos Yáñez. Estadística Espacial con R. https://rubenfcasal.github.io/estadistica_espacial/index.html 2023.Paula Moraga. Spatial Statistics for Data Science: Theory and Practice with R. https://www.paulamoraga.com/book-spatial/index.html Mariano Córdoba, Pablo Paccioretti, Franca Giannini Kurina, Cecilia Bruno, Mónica Balzarini. Guía para el análisis de datos espaciales. Aplicaciones en agricultura. https://www.agro.unc.edu.ar/~estadisticaaplicada/GpADEAA/ --- class: inverse middle